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ABSTRACT 



The problem of estimating the position of an extreme 
point of an unknown function of several independent variables 
is examined for the case where the dependent variable is known 
to be bounded. The classical method of least sequres is formulated 
as a quadratic programming problem to be solved numerically on a 
digital computer, where the coefficients of the fitted equation 
are determined subject to restrictions on both the independent 
variables and the dependent variable. 

Several two dimensional models were examined using synthetic 
experimental design techniques. The results, though not conclusive, 
indicate that the method of bounded least squares can be a useful 
computational tool in some two dimensional problems. It remains 
to be shown whether the algorithm is useful in problems involving 
more than two independent variables. 
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INTRODUCTION 



The method of least squares to fit a curve or surface to a set 
of data points is well known in the literature, as are techniques 
for determining the extreme points of such a surface. In par- 
ticular, non-linear methods are available £loJ to examine the 
surface in a particular region defined by bounds on the independent 
variables. 

While considering the special problem of approximating, by 
a quadratic equation, a complex probability distribution, it was 
proposed to make use of optimizing techniques which would give a 
good estimate of the particular combination of independent variables 
producing maximum probability. Since the equation was only an 
approximation, it was expected that errors would occur in the 
estimation and that further accuracy would be obtained by testing 
in the physical system. 

At this point it occured to the author that a large class of 
problems of this type must exist; i.e., certain independent variables 
are known to be bounded, and it is desired to find an extreme point 
of some function of these variables, which is itself bounded. 

A considerable body of literature on the theoretical nature 
of non-linear programming exists to consider the bounded independent 
variables M H , and it seemed reasonable to extend these 
methods to treat the bounded dependent variable in such a manner that 
the extreme point of the true Response Surface could be more closely 
estimated. 
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STATEMENT OF THE PROBLEM 



Consider a physical process or system under investigation 
which has at least one characteristic called the Response, such 
as cost, yield, power, etc. Assume that this Response is some 
unknown function, called the Response Function, of a set of 
independent variables and that these variables can each be quantita- 
tively fixed at a given value, or level. The observed response 
resulting from a combination of the variables, however, is the 
sum of the true Response and a random error which is distributed 
about a mean of zero with constant variance (no matter what the 
values of the system variables). 

Geometrically, the Response Function may be represented by a 
surface in (n+1)- space called the Response Surface. The ex- 
perimenter has systematically proceeded to a particular region, 

R, in n- space (defined by the n independent variables) which is 
bounded by the planes x^ ~ d^ and x^ ^ g^ (i- 1, 2, . . ,n) . He then 
wishes to find the levels of the independent variables for which 
the Response is extreme, restricting the independent variables to 
lie within R # 

To do this he must make the assumption that the Response 
Function may be adequately represented by some analytic function. 

In our case we shall assume that a quadratic function of the form 




(i) 



2 



will adequately represent the Response Function, at least in the 
vicinity of an extreme point. 

If we also assume that the experimenter has some prior know- 
ledge that the Response Function is bounded, then he should be 
able to make use of this knowledge in any procedure which estimates 
the coefficients of the quadratic equation used to locate the extreme 
point. 
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CLASSICAL METHOD OF SOLUTION 



In a classical solution to a problem of this type, experiments 
would be performed at certain selected points in the region, R, 
and the method of least squares £9] would then be employed to 
determine the coefficients of the equation to be fitted. 

Once the quadratic has been estimated, the set of partial dif- 
ferential equations derived by differentiating the function, Y, 
with respect to each x^, can be solved for the location of the 
extreme points. In order to remain within the region defined by 
the bounded independent variables, nonlinear programming methods 
may be applied to determine the position of the bounded extrema £lo] . 

We notice that no restriction is placed on the Response, Y. 

The optimizing technique has given the set {x^*; which cor- 

responds to an extreme point, Y -H(x^ ,...,x^ ), of the fitted 

<jS* 

equation only. It is quite possible that Y does not lie within the 
a priori bounds on the true Response and, in practice, further ex- 
perimentation could be performed at and near |x^ , ...,x J ** n or< * er 

to confirm the existence of the extreme point and to more exactly 
determine the Response at that point. If the research worker has 
the time and resources available to conduct the large number of 
experiments necessary to fit the quadratic by the method of least 
squares and then to converge on the true position of extreme response 
from a predicted position, such a method may be entirely adequate. 
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It was felt, however, that such a procedure was not as efficient 
in predicting an extreme point as should be possible when a priori 
bounds on the Response are known. Accordingly, the classical method 
was examined to determine where a modification to the theory could 
best be made. 
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THE METHOD OF BOUNDED LEAST SQUARES 



Though further experimentation is to be performed in the 
vicinity of the predicted extreme response, the problem of optimization 
is concerned not only with predicting response, but with predicting 

'ft's? 'f c 'f c 

the Y -H(x^ ) with sufficient precision that the least 

number of experiments are necessary to converge to an extreme point. 

In the classical solution, two options are available. First, 
the bounds on the Response may be ignored and the position of maximum 
response obtained from the fitted equation may be taken as the best 
estimate of the values of the independent variables which will maximize 
(or minimize) the observed response. Secondly, the hyperplanes 
representing the upper and/or lower bounds on the Response may be 
intersected with the Response Surface obtained by the estimate of 
the quadratic equation. Only those portions of the fitted response 
surface lying within the bounded region would be searched for an 
extreme point (see figure 1), This method was not considered to 
yield an acceptable solution since the entire range of values of 
the independent variables corresponding to the portion of the fitted 
response surface outside the bounded region would give an extreme, 
restricted response. 

Situations can exist in which the fitted response surface 
varies greatly from the true Response Surface, both in shape and 
range of value over the surface. Thus, it was hypothesized that 
the extreme point of the fitted response surface could be moved 
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closer in n- space to the extreme point of the true Response Surface 
if the quadratic equation used to represent the true Surface could 
be estimated subject to the restriction that no point on the fitted 
surface could fall outside the a priori bounds on the true Response 
(see figure 2). 




TRUNCATED LEAST SQUARES FIT 



Y 




BOUNDED LEAST SQUARES FIT 



Recall that equation (1) defined the fitted, or calculated 
response. If Y is defined as the observed response, then 
(k 5 l,...,m) are the sets of response values obtained when m dif- 
ferent combinations of values of the x. are 1) set as levels in 

l 

I 

experiments where Y is observed and 2) used to calculate Y from 
the fitted equation. Then, 
m , 

G * X (Y, - Y,) 2 (2) 

k= 1 
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is the sum of squares of the residuals. The sura, G, is sometimes 
referred to simply as the Residual. 

Thus, the classical least squares method of estimating the 
coefficients may be modified as follows: 



Minimize G with respect to all coefficients of the x^ (3) 
subject to &£ * % i * (i s l,...,n) (4) 

and e ^ ^ f (k- 1, . . . ,m) (5) 



where n is the number of variables and m is the number of experiments 
(data points). In Appendix A it is shown that an equivalent form 
of (3) is to Minimize F^ with respect to all A^ , where 



F * z 



i-1 



CoA. 

i i 



j-i 



CK i A i 



( 6 ) 



and that (5) may be restated as 
Pj * 0 (j-2 2mtl) 



(7) 



Non-linear programming methods designed to accept problem state- 
ments in the above form, where the constraints are linear in the 
unknowns, are classed as quadratic programming. The particular 
algorithm used in this paper is due to Klingman^ £gj, and was selected 
for its general application to the case where non-linear constraints 
are encountered in addition to bounded Response, as well as its avail- 
ability in finished program form. 



The program originally published in Klingman's thesis was 
modified after private communications with W. R. Klingman and Dr. 
D. M. Himmelblau, University of Texas, Austin, Texas. 
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EXPERIMENTAL PROCEDURE 



Almost any convenient and reasonable method may be used to 
select the points at which the experiments are to be performed. 

In keeping with the image of an economy-minded experimenter, how- 
ever, it was decided to make use of the techniques of composite 
design due to Box £lj £4j. As shown in Table I, a considerable 
reduction in the number of points to investigate is possible when 
the 2° factorial design is augmented rather than using the 3° 
factorial design usually employed when quadratic effects are being 
observed. 



TABLE I 

COMPARISON OF COMPOSITE DESIGN 
AND FACTORIAL DESIGNS 



no. of 


no, of exi 


jeriments required 


variables 


composite 

design 


3° design 


2 n design 


2 


9 


9 


4 


3 


15 


27 


8 


4 


25 


81 


16 


5 


43 


243 


32 


6 


77 


729 


64 


7 


143 


2187 


128 


8 


273 


6561 


256 



In order to simulate actual experimentation, a hypothetical 
physical system was represented by an analytic equation. For each 
combination of independent variables selected by the design program 
in Appendix B, the Response was evaluated and a random error dis- 
tributed with zero mean and constant variance was added to it to 
obtain (k s l,...,m). Using the data points thus generated, the 
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coefficients of the quadratic equation were estimated by 1) the 
classical method of least squares [6] M and 2) the method of 
bounded least squares (Appendix B) . This procedure was repeated 
until a total of ten estimates of the equation were obtained for 
a constant response error variance, for each of four analytic 
equations tested, the independent variables being set at the same 
levels for both methods. Each of the fitted equations was then 
solved for an extreme point in the design region. The distance, 
in n-space, between the true extreme point and the extreme point of 
each fitted equation was determined. Finally, for each case, consisting 
of ten runs for a particular equation, the distances were averaged 
to obtain the average extreme point error, r. 
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RESULTS 



For each case, the distance in n-space between the true and 
estimated position of extreme Response was found using the relation 




and then averaging over the number of runs to obtain r. This 
procedure was followed for the results of the classical least 
squares fit and the bounded least squares fit for each case. 

Numerical Results 

CASE ONE j Y - 4 x ][ x 2 - - x^ 

The response error variance was set at .25 and the region, R, was 
defined by^( x ^ • 0-x^2, O^x^- 2^ . A maximum is located at (1,1) 
and an a priori upper bound on the Response was assumed to be 1.0 . 
Using the method of bounded least squares, the estimated maximum was 
located 12.7 percent closer to the true maximum than by the classical 
method. 

CASE TWO : Y s -x 4 4 2x 3 4 2x 2 - 2x + 10 
The response error variance was set at 1.0 and the region, R, was 
defined by ^x; ~2 ^ x - 3^ . An a priori upper bound on the Response 
was assumed to be 15.0 and two maxima are in R; a local maximum at 
x =*-1.0 and a global maximum at x - 2.0. Using the bounded least 
squares method, the estimate of the extreme point was only 1.2 per- 
cent closer to the true maximum than by using the classical method. 
This difference was found to be completely submerged in the effect 



i 



of response error variance upon the location of the maximum 



TABLE II 

COMPARISON OF MEAN EXTREME 
POINT ERROR FOR CASE ONE 





Classical Method 


Bounded Method 


RUN 


kk 

x i 


kk 

X 2 


r 


kk 

x i 


** 

X 2 


r 


10 


0.7007 


0.6103 


0.4914 


0.7188 


0.6403 


0.4566 


9 


.7614 


.6518 


.4219 


.5721 


.5986 


.5867 


8 


.7607 


.6748 


.4037 


.6647 


.6228 


.5047 


7 


.6593 


.5934 


.5304 


.5177 


.5758 


.6423 


6 


.6403 


.6133 


.5282 


.6961 


.6568 


.4585 


5 


. 7932 


.6993 


.3647 


.6312 


.5483 


.5832 


4 


.7650 


.6494 


.4219 


.6181 


.6277 


.5334 


3 


.6004 


.5493 


.6022 


. 7293 


.6308 


.4578 


2 


.0000 


.0000 


1.4142 


.5742 


.5389 


.6277 


1 


.2386 


.4357 


.9476 


.7303 


.5787 


.4970 


r 


.6126 


.5348 



CASE THREE : Y - 



(x r + l) 2 + x 2 2 + 10“ 8 



The response error variance was set at .0001 and R was defined 
by s ~4 - x l * 4 S 4 ° x 2 - 4^ . The Response Surface 

is symmetrically oriented about a maximum located at (~1,0), and 
an a priori upper bound of 1.0 was placed on the Response. Both 
methods of fitting the experimental points resulted in identical 
estimates of the maximum at (0,0) for all runs. It should be noted 
that the Response Surface is very shallow outside that area enclosed 
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by a cylinder of radius 1.0 which is centered at (-1,0) 



TABLE III 



COMPARISON OF MEAN EXTREME 
POINT ERROR FOR CASE TWO 





Classical Method 


Bounded Method 


RUN 


** 

X 1 


r 


** 

x i 


r 


10 


0.5669 


1.4331 


0.5706 


1.4294 


9 


.5603 


1.4397 


.5670 


1.4330 


8 


.5488 


1.4512 


.5853 


1.4147 


7 


.5603 


1.4397 


.5623 


1.4377 


6 


.5636 


1.4364 


.5897 


1.4103 


5 


.5431 


1.4569 


.5670 


1.4330 


4 


.5540 


1.4460 


.5735 


1.4365 


3 


.5400 


1.4600 


.5809 


1.4191 


2 


.5702 


1.4298 


.5537 


1.4463 


1 


.5636 


1.4364 


.5980 


1.4020 


r 


1.4429 


1.4251 



CASE FOUR : Y - Xj 2 - - 3x^2 

The response error variance was set at .01 and an a priori lower 
bound of -1.69 was placed upon the Response. The region, R, was 
defined by £(x^,X 2 ): - l-x^-3, -1-X2~2| . There is a saddle point 
at (0,0), a minimum at (2.25,1.5), and the global minimum for R is 
located at (-1,-1). Both methods correctly estimated the position 
of the minimum for all runs. 
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Computational Results 



Program DESIGN required about one minute to compile and 
execute for two variables, including printout and punch. Program 
FIT required approximately three minutes to compile and four seconds 
to execute for each run. The minimum fractional change in both the 
coefficients being estimated and the value of the Residual being 
tested was set at 10 In his thesis, Klingman gives a compile 
time of about one and one-half minutes for the quadratic program 
(subroutines DIRECT, PATSER, EXPSER, GRAD, SUMS and OUTPUT in program 
FIT of this paper). 

The program due to Dillon was used to fit the quadratic by 
the classical method. That program consists of the elements of 
programs DESIGN and FIT with the equations from Krase and Cyl-Champlin 
L»3 substituted for the subroutine package which comprises the 
quadratic program portion of FIT. The program required one and 
one-half minutes to compile and less than one second to execute for 
each run. Thus, the classical method of solution required about 
one-half the total computer time required by the bounded least 
squares method. As mentioned before, no attempt was made to reduce 
either compile or execution time, for either program, to a minimum. 
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CONCLUSIONS 



The algorithm presented here was expected to provide a closer 
estimate, in terms of the physical values of the independent variables, 
of the location of a point of extreme Response. It was assumed that 
the method would be of greatest value in cases where the Response 
Surface 1) is non- symmetrical about an extreme point in the interior 
of the region being investigated or 2) has more than one extreme point 
in the same region. 

In no case tested did the algorithm provide worse estimates of 
the true extreme point than did the classical method. For cases 
three and four, where the Surfaces were symmetric, and oriented with 
an extreme point at the boundary, no particular advantage accrued 
from using the method of bounded least squares. In case two the 
results were inconclusive, and case one did not produce the desired 
improvement over the classical method, even though a distinct im- 
provement was observed. 

For one and two dimensional problems, it is felt that the 
advantage of the system is outweighted by the disadvantage of 
observing the Response at only a small number of points. Any 
method involving least squares requires enough points to adequately 
represent the Surface being approximated. In the bounded least 
squares method, where the advantage lies in estimating an extreme 
point of a complex surface, composite design point selection does 
not choose enough points to adequately define the surface. Further 
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test® should he made to determine the minimum, number of points 
necessary to permit the method of bounded least squares to enjoy 
a distinct advantage over classical methods, 

No suitable problems of dimensionality three or higher were 
found which permitted rigorous testing of this algorithm. Since 
the number of experimental points required for these problems 
becomes quite large (see Table I) p a significant advantage may 
be realised by the method of bounded least squares, since additional 
information is included in the solution to the problem at no ex- 
perimental cost. The cost of further testing, as opposed to doubling 
the computer time, must be considered for each system investigated, 
before the decision is ma.de to select any method which reduces 
the number of experiments. 
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APPENDIX A 



FORMULATION OF THE METHOD OF BOUNDED LEAST 
SQUARES AS A QUADRATIC PROGRAMING PROBLEM 



All terms used In this Appendix have been previously defined 
except the following: 

£ : magnitude of the bounds on the coded values of all the 
independent variables, X^. 

K : constant due to orthogonality requirement which permits 
independent evaluation of each coefficient. 

: coefficients of the orthogonal quadratic 
equation, for coded 

(NOTE ; double letter symbols enclosed by ( ) refer 

to names of single variables.) 

The form of the quadratic equation with which we desire to fit 
the experimental data is 




Y 




z 

1=1 



B l x i + 



£ 

1=1 



(BB) ii (X i 



K) 



+ £ I (BB) X X 
I s 1 j>i 1J 1 J 



(A- 1) 



so that the Siam of squares of the residuals may be writted as 



n 9 

-z 



k= 1 

m 



V 



— r - > 

s 2- <\> 
k=l R 



Y, 



k* 1 



±< 



k- 1 



k Y k 



(A- 2) 



Simplifying each term, 





® constant 



(A- 3) 



19 



and 



m. 2 m n.. 

B » * li li 0 ^ 



m n , 

* Z 21 (BB) (X - K) 

k=l i-1 11 1 



((BB) lj X ik X jk ) 2 



m n jq^ 

+ III 

k- 1 i»l j>i 

where all other terms can be shown to be zero. Then, using the 
relations developed by Krase [V] for classical least squares 
theory using orthogonal designs, 

1 m 



(A-4) 



* — X x 



m k“ 1 



ik 



for any i=l,2,...,n 



“ — (2 n - li 2 ) - JF 

m 

<£ 2 = - 2 n ) 



(A- 5a) 



(A- 5b) 



(A- 5c) 



m 



kn 



m 



(X 



ik 



K) 2 - 2<$^ for all i“l,...,n 



21 x 2 x ik 2 s 2° foff a11 iVj 

k-1 1 J 



m 



X (x lk 2 



K) - 0 for all i*l,...,n 



k*l 

it is easily shown that 



m 



X Y. 

k~l 



2 - „ 2 
k " mB o 



+ mK ^ B ' 
1*1 



(A-5d) 



(A-5e) 



(A-5f) 
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(continued) + 2 (BB) 



i=l 



n n 



ii 



+ 2 n 2. 2. ( BB > 

i a l j>i 



ij 



(A-6) 



The last term of (A- 2) is a straightforward sum of multiplica- 
tions. When terms are combined. 



g=z j b . 2 . 2 <5 4 J*" 

s fh. 1 m 

+ 2 n £ £ (BB) .. 2 - 2 21 Y.B^ 



(BB) 



ii 



i-1 j-i 



' Ii Bl { 2 Ii <x ^i 
' Sr <BB> “ 1 2 Si Y ^ <x ^ 2 • K> } 



n n 



i=-l j-i 



m 



(BB)..(2 > Y,X.,X.,) 



Then simplifying 



ij 



n 



^ k ik jk y 



(A- 7) 



n n 



G - Z + Z B 
s o 



O+I Z i B i +1 I (Z2) i1< BB >i1 

i= 1 11 i»l 1»1 1J 1J 



2 2 ^ ^ 

+ V. + Vi + 



£ ij <TO lJ (BB) lj : 



(a- 8) 



where the substitutions are obvious. 

The Response is bounded by an upper and lower restriction 
which can respectively approach plus and minus infinity; 

( k— 1,2,..., m) 

Then, (A- 9) can be restated as 



** Y k* f 



Y, - e - 0 



f - Y, - 0 



( k* 1,2, ... ; m) 



(A- 10) 



21 



or, using (A~l) and the trans formation 



\ = X A !<>•■*> lk (A- 11) 

a change of variables may be made so that 
F - 0 (i-2,3, . , . s ro*l) 

1 (A- 12) 

Fj - 0 (j=m*2, . . . ,2ra* 1) 

Similarly, (A“8) may be transformed to 

(CK)^ 2 (A- 13) 

i-1 

Equivalent quantities for the case where the number of 
variables, n, is equal to two are listed in Table IV. The 
equations (A- 12) and (A~13), together with the bounds on the 
Xp constitute the form of the quadratic programming problem 
solved in Appendix B. 



kn 



G + 2^ 4 

B 



i“l 



c i A i + 



TABLE IV 

EQUIVALENT QUANTITIES IN THE QUADRATIC PROGRAMMING PROBLEM 



FOR TWO VARIABLES 



Unknown 


Objective 




Constraint 


Constants 


Coef f ic tents 


Constants 




(k = 1,2, 


. . . ,m) 


B o 


A 1 


Z o 


C 1 


rr H 

et n fp 


1.0 


“l,k 


B 1 


a 2 


Z 1 


C 2 


o m 
m ^ 

!% m n 


X l,k* ‘ 


“2,k 
















b 2 


A 3 


Z 2 


C 3 


H* 3 © 

Q H* 3 

^ H ® 
0 

M H H* 
* O 

ft H 


X 2,k 


“3,k 


BB U 


A 4 


zz n 


C 4 


x i% - K 


n 4,k 


n o 


BB 12 


A 5 


zz 12 


C 5 


o 

n 


X l,k X 2,k 


^S.k 


BB 22 


A 6 


zz 22 


C 6 




X 2% - K 


“e.k 
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APPENDIX B 



The FORTRAN code for selecting the experimental points 

(program DESIGN) and for finding the coefficients fo the quadratic 

equation by the method of bounded least squares (program FIT) is 
2 

described here. The production runs were made on a CDC 1604 
digital computer under COOP Monitor control. The programs allow 
considerable generality as written, and no attempt was made to 
"clean up" either one for the specific use of a bounded least squares 
fit. 

A general flow diagram is not considered essential to an 

understanding of either program DESIGN or FIT. Subroutines DIRECT, 

PATSER, EXPSER and SUMS comprise the Multiple Gradient Summation 

algorithm [8j used to solve the quadratic programming problem. 

Detailed flow charts for these routines are available from other 
3 

sources. 

As presently programmed, a maximum of seven variables can 
be considered. On a 32K machine additional variables can be 
included if only one bound (upper or lower) is placed on the 
Response. Subroutines SUMS and GRAD must then be modified 
appropriately. 

^The algorithms given in j^6^| and [8] were modified and combined 
for this use. 

3 

Dr. D. M. Himmelblau, The University of Texas, Chemical 
Engineering Department, Austin, Texas, 
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1 . FORTRAN Nomenclature . 



Only the essential symbols used in the FORTRAN programs of 
this appendix are listed. 



A (I) 


- The coefficients to be found. 


AK 


- Orthogonality constant (K in Appendix A). 


ALN(I) 


- The best set of coefficients found. 


ALO(I) 


- Previous base value of A(I). 


BZRO,B(I), 

BB(U) 


- Coded coefficients of the fitted response function. 


BO 


- Coefficient (see REALB, REALC) . 


C(II) 


- Renumbered constants. 


CINE(I) 


" Variable used to determine first negative 
constraint crossed. 


CK(U) 


- Renumbered constants. 


DEL(I) 


- Incremental move for each A (I), 


BELT 


- Bound on the coded variables (<£ in App. A). 


EP(I) 


Incremental change in A (I). 


EV 


- Estimate of error variance of the fitted 
response values at the design points. 


EX(I,J) 


- Physical value of independent variable I 
in experiment J. 


F(I) 


- The objective function (Residual) and the 
constraints. 


FCHG 


- Minimum allowable fractional change in F(l) 
before a search failure occurs. 


FN 


- Value of F(l) tested in types I and II 
exploratory searches. 


FNN 


- Value of F(l) calculated in subroutine SUMS. 
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FNEG(I) 


The negative constraints when a non- feasible 
initial guess of A(X) is made. 


G 


- Sum of squares of the residuals. 


XACCEL 


- Accelerating factor in the adaptive move. 


INK 


Frequency for printing values during a successful 
pattern move sequence. 


ILA-TE(X) 


Latest constraint contacted during an adaptive move. 


ITED 


Search failure counter. 


ITER 


- Pattern move failure counter. 


K 


Used only in DIRECT and associated subroutines 
as the number of unknowns (coefficients) to find. 


KCON 


Total number of constraints contacted in an 
adaptive move. 


KDAPT 


- Counter for the number of adaptive moves. 


KFREQ 


Pattern move/ exploratory search failure printout 
frequency. 


KFREQS 


- Successful pattern move printout frequency. 


KN 


Number of coefficients of the fitted response 
function (also number of each type of coefficients 
C and CK) . 


LX»LY 


If equal to 2 S the number of search and pattern 
move failures/successes will be printed. 

If equal to l 9 no printing will be done. 


M 


= Number of experiments. 


MK 


- Counter for exploratory search failures. 


MM 


Counter for successive exploratory search failures. 
Used to reduce step size. 


MQ 


- Maximum number of allowable successful pattern moves 


N 


- Number of independent variables (only for main 
programs DESIGN and FIT and subroutine REAL)* 
maximum of seven. 
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Number of constraints plus one (only for sub- 
routines DIRECT, PATSER, EXPSER, SUMS, GRAD and 
OUTPUT) - equals 2M ■+ 1. 


NC 


Identifying case number. 


NOMEX 


Counter for search failures. 


NOPTN 


Counter for successful pattern moves. 


NUMEX 


Counter for number of exploratory searches. 


NX 


Number of extra points at which to evaluate the 
fitted response function. 


PARTIAL(I,J) ■ 


- Partial derivatives of F(I) with respect to A(J). 


PER(I) 


DEL (I) equals FER(I)*ALO(I) . 


PERC 


Initial fractional incremental move for each A(I) 


QO,Q(I), 

QQ(J) 


Constants as defined in Appendix A, 


REALB(I) , 
REALC (I, J) - 


Coefficients of the fitted response function 
when physical values of the independent variables 
are used. (see BO also) 


RL, RU 


Lower and upper bound on the Response. 


SLOPE (I) 


Factor used to convert to coded values of the 
independent variables. 


SUD(I) 


Number of successful successive search moves 
for each variable. 


TEST 


If PER(I)^TEST, the problem has converged to 
a final answer. 


VECTOR(I) 


Sum of gradients of the contacted constraints 
in an adaptive move. 


WAY 


Set equal to -1.0 to indicate that a minimum 
of the objective function is desired. 


X(I,L) 


Coded value of the independent variable I in 
experiment L, 
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XAV(I) 


- Constant used to convert to coded values of the 
Independent variables. 


XMAX(I) 


- Upper bomd on the EXCX^J) for all J. 


XHIN(I) 


- Lover bound on the EX(I y J) for all J, 


XX{I 9 L) 


* Constants as defined in Appendix A. 


Y(L) 


“ Observed response (except in subroutine PATSER). 
Temporary storage in subroutine PATSER. 


YD(L) 


- Residuals for each experiment. 


YY(L) 


- Computed value of the response function. 


Z0 9 Z(I), 

ZZ(J) 


Constants as defined in Appendix A. 
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2. Input for Program DESIGN . 



Card One : 


FORMAT (212) N, NC 


Card Two : 


FORMAT (6E12.8) . . » .XMIN(I), XMAX(I) 
punched in pairs for each variable, I - 1,N. 
Maximum of three pairs per card. 



. Input for Program FIT . 



Card One : 


FORMAT (215,E12.8) . .N, M, DELT 


Card Two : 


FORMAT (6E12.7) . . . .SLOPE(I), XAV(I) 
Punched in pairs, I - 1,N. May require up 
to three cards for this data. 


Next Cards; 
(total of 
M cards) 


FORMAT (8E10.3) . . . .X(I,L) 
up to eight values of the independent var- 
iables for the Lth experiment. 


Next Cards: 


FORMAT (8E10.6) . . . . Y(L) 

observed response for up to eight experiments 
per card. 


Next Card : 


FORMAT (215,2E12.8) . .NC, NX, RU, RL 


Next Card : 


FORMAT (80 character Hollerith title) 


Next Card : 


FORMAT (15,F14.8) . . .MQ, PERC 


Next Card : 


FORMAT (8E10.6) . . . .Initial guesses 
for values of the A(I). 


Next Card ; 


FORMAT (E10.6,215). . .TEST, LX, LY 


Next Card : 


FORMAT (E 10. 6, 215). . .FCHG, KFREQS, KFREQ 


Last Cards: 


FORMAT (8E10.6) . . . .EX(I,NX) 

Extra points to evaluate. Punch I independent 
variable values per card; NX cards. 
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PROGRAM DESIGN 

DIMENSION X(8, 273) , EX(8,273) 
COMMON/BLOCK2/SLOPEC 8 ) ♦ XAV(8) 
COMMON/BLOCK3/XMIN (8 ) » XMAX(8) 
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1116 FORMAT ( 4H 1 N=,I2,4H, M= , I 4 , 7H , DELT A = , E 1 4 . 8/ 93H0THE EXPERIMENTAL VA 
1 LUES OF THE INDEPENDENT VARIABLES WITH THE CORRESPONDING CODED VAL 
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C**********EQUATE Z AND 0 TO SINGLE SUBSCRIPT CONSTANTS********** 
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885 FORMAT (215, 2E12.8) 



C**********EQUATE coefficients from quadratic program to standard form 
c for use in subr. real in order to use real VALUES OF X ( I ) 
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PRINT 100* (NC* N, SIGMA, M, RSQ, KN, EV, DELT, G, N2» ( 

1 L » X ( 1 » L ) » X ( 2 » L ) » X ( 3 » L ) » X(4,L)» X(5,L), X(6»L), X(7,L), X(8,L), 
2 Y( L ) * YY ( L ) , YD ( L ) , L= 1 * M ) ) 
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90 CONTINUE 

READ 2 » MO, PERC 
2 FORMAT ( I 5 » F14. 8 ) 
PRINT 105, PERC 
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50 FORMAT (28H PROBLEM FAILED TO CONVERGE. ) 

PRINT 82 >NUMEX > NOPTN 

82 FORMAT (30H NUMBER EXPLORATORY SEARCHES. = 15// 39H NUMBER OF SUC 

1CESSFULL PATTERN MOVES = 15/) 

CALL OUTPUT (A, K) 
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80 CONTINUE 

DO 85 L=1 »M $ LP1 = L + 1 $ LP1M = LP1 

85 F(LP1M) = -F(LP1) + RDIFF 
FNN = F ( 1 ) 

ISTOP = 2 



CO 

•> 

vO 

e. 

vO 

*•— < 
u_ 
I 



















z 
























































o 












3 












• — • 










4 




•> 














f\J 




u 












M 












3 














< 














»— 1 
















< 












W. 














w 












‘ — • 


















Ll 












O 


























*— 1 


CO 




UJ 






*— t 










co 








LU 














o 














•> 




h- 




— 








•s 












Z 














< 


vO 










Ll 


CO 




< 






II 






• — • 




co 






LU 


Ll 














(X 


•* 












— 




-J 




— 








f— H 




co 






z 
















o 


vO 










II 






t—4 




Ll 








CO 




r. 






1— < 


+ 
















•« 












CM 










z 






•> 




CO 






















LU 


in 




rH 


I — t 


rH 








UJ 




ll 


o 




»— H 


o 




CO 




o 


»-H 


Z 










l-H 




z 




UJ 








i— H 


1 


UJ 


z 


UJ 




u 


UJ 


1 


CO LU 


. — . 


LU 


• 


II 


Z LU 


LU 




3 






l-H 




3 


+ 


II 


II 


1 




3 


*— 1 


3 






3 


3 




3 


LU 


3 O 


3 U_ 


3 


3 r- 


3 


n 




h- 


• — » 


Z 








3 


3 


Z 


u 


Z 


►— < 




Z 






Z 


z 


Z 








Z 


z 




Z 


z 




3 




*—• 


3 


UJ 


a 






*— * 




•— • 


— * 


LU 


►— i 


II 


O 


*— i 


• — « 


♦— i 


II 


CM 


II 


i— < 


* — « 


O 


*— « 


KD C£ 




O 


Ll 


h~ 




z 


o 


o 


W- 


h- 


w. 


h- 


LU 


h~ 


h- 




h- 


1— 




►— 




CO 




h~ 


1— 


h- 


H- 


X 3 




(X 




Z 


II 


►—I 


h- 


UJ 




z 




z 


Z 


< 


z 


UJ 




z 




z 


z 




z 


z 


z 




z 


U h- 


a 


co 


Ll 


o 




V3 


CO 


z u. 


o 


Ll 


o 


l-H 


3 O 


ZOOlL 


o z 


o z 


o 


O O O U- 3 


z 


3 


i — i 


<3 


3 


»— 1 


»— 1 


Ll 


M 


V3 


*— i 


o u 


*— * 


3 


*— * 


o c 


• — t 


u 


Ll 


o u_ 


uuou 


♦-» ££ LU 


lO 




in 












CJ 




CO 






vO 






o 












CM 


rH 




CO 


r- 




















«-H 
















CO 




CO 








CO 


CO 




CO 












DIMENSION A(A5) 



cl 

3 

CL 

</) 

rsj 



in 



< 

li- 



lt \ 

s 

in 



— o 

»-H 

H-l • 

— O 
z <nj 



vO 
































in 




< 






co 
































<3- 




•*— 






































— C\J 




Ll 






r- 
































O CL 










00 
































3 3 




z 






CM 




























* 




m 3 




o 






w 
































* o 




►— < 






3 




































h- 


n 


•— t 


< 
































in »~h 




z 






•— i 
































<t’CL 




3 


m 


3 


h- 
































— 3 




Ll 


*—* 


UJ 


CL 
































Z 13 








o 


< 


























\ 






3 O 




UJ 






CL 






























< * 




X 




— 




























s: 






* CNJ 




t— 


i— i 


»— I 




























o 






— < 








•*— 


in 


























z 






m cl 




Ll 


< 


CL 


<3* 


























13 






<3“ < 




O 


>— 


3 


— 




















<\i 






O 






— a 






Ll 


CL 


X 














— 






* 






X 






o ~ 




13 






U 














►— 4 






* 












3 ♦ 




3 


Q 


















* 






— 












< -I 




3 


UJ 


















r-H 












3 






» < 




< 


l— 


*— • 


in 




— 


* 








CL 






3 






#> 






— CL 




> 


m 




<3- 














3 






* 






►— i 






m < 






UJ 


< 






X 
















*— • 












<3- CL 




UJ 


h- 




U 




u 










3 


















•*— *—* 




m 






* 




* 










< 


















3 * 




< 


z 








• 










*— < 












3 


— 




13 s : 




CD 


o 


M 


co 




OvJ 






3 




t— 






3 


— 




< 






O ST 






*— I 


**«-' 


<* 










•» 




CL 






< 






t-H 






•s •> 




1— 


1— 


z 


t—i 




+ 










< 






*— • 


s: 




1— 


* 




— z 




Z 


u 


3 












>— 




CL 






i — 


o 




CL 


< 




m u_ 




3 


z 


< 


vO 










X 




1 






CL 


z 




< 






<r * 




CL 


3 


* 


CO 




•— i 






X 










< 


UJ 




CL 






— o 




CL 


Ll 




— ' 














II 






CL 


Q 






h- 




CL 3 




3 




— 


X 




u 






II 
















II 


3 




13 O 


Z 


U 


UJ 


i— » 


X 














— - 






+ 








a 




CL U- 


Ll 




X 




X 




II 






















<— ■> 


i— 


in 


X 


* 


3 


f— 


O 


On 




















s: 


Ul 






3 




h- 


O 


X 




3 


X 


Z 


— - 








s: 


s: 


Z 




X C 


t— 


x 


3 


O 


— 


X 


3 


1 — 


X 


< 


u 


X 


HH 


s: 




r— 1 


+ 


?— i 


•> 


o 


* Z 


CL 


•s 


* 




< 


U 


O 




<3 




o 










CL 


»— i 


a 


r-H 


• 


r*H 13 


O 


►— i 


13 




o 


Ll 


X 


co 


X • 


3 


f— H 


r~t 


«— i 




3 


CL 


3 




o 


Q 


in 






Z 


z 


3 




<3- 




O' « 


CD 


II 




ii 


+ 




3 


— 3 


II 




II 


3 


II 


13 




o 


CD 


O 


in 


X 


»— 1 



o 

u 



< 

oo h- 
CL 

o < 

O CL 



II 

oo 

»— i 

O CL 
O -J 



JD II II 3 II -J 3 I— •— X 

< z -> z 3 < z z 3 in z 
*--- s: s: oc ozo 

hh<r o^oi-o^i-hD cl uj s: 

Z ZZZ QCZI-Q costs: 
UJOUJOLUOCOUJZ 3 O 

ooouooduo^iu m o u 






CO II 



00 



CL Z 

< o 

CL U 



h- < 

Z Z tZMH 

MQ1 O M CL 

CL O ^ CL O d 

CL Ll UJ CL O 



<3* CO <t 

oo 



o c\j <f co 

O' 



52 



FORMAT ( I3,5E23. 10) 

FORMAT ( 3H AI,7X»1AH PREVIOUS BASE 10X.13H CURRENT BASE 1 IX * 1 OH TEST 
BASE1 2X ♦ 15H PERCENT CHANGE12X »6H DELTA //) 

END 



APPENDIX C 



SAMPLE COMPUTER PRINTOUT FOR CASE ONE 
N= 2 *M= 9 ,DELTA=1.00000000E 00 , 

THE EXPERIMENTAL VALUES OF THE INDEPENDENT VARIABLES WITH THE 
CORRESPONDING CODED VALUES ARE 

FOR INDEPENDENT VARIABLE NUMBER 1 THE EQUATION FOR THE 
DESIGN POINTS IS EX=AX+B, WHERE 

EXM I N= 0 AND EXMAX= 2. OOOOOOOE 00 

A= 1 • OOOOOOOE 00 AND B= l.OOOOOOOE 00 



EXPERIMENT 


- X CODED 


EX REAL 
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.9 
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l.OOOOOOOE 
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FOR INDEPENDENT 


VARIABLE NUMBER 


2 THE EQUATION 


FOR THE 


DESIGN POINTS IS 


EX=AX+B ♦ WHERE 







EXM I N= 


0 AND EXMAX 


= 2. 


OOOOOOOE 00 




A= l.OOOOOOOE 


00 AND B 


= 1. 


OOOOOOOE 00 
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00 
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